﻿#include "ComputeArea.h"
#include "FiniteElement.h"
#include "Quadrature.h"

#include <vtkm/cont/Algorithm.h>
#include <vtkm/worklet/WorkletMapTopology.h>

#include "Mesh.h"

struct CalAreaWorklet : public vtkm::worklet::WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn, FieldInPoint, FieldOutCell);
  using ExecutionSignature = void(_1, _2, _3);

  template<typename Fty, typename... Args>
  VTKM_EXEC void OnCellOperator(vtkm::CellShapeTagGeneric tag, Fty&& f, Args&&... args) const
  {
    switch (tag.Id)
    {
      case vtkm::CELL_SHAPE_TRIANGLE:
      {
        using CellShapeTag = vtkm::CellShapeTagTriangle;
        f(CellShapeTag{}, std::forward<Args>(args)...);
      }
      break;
      case vtkm::CELL_SHAPE_QUAD:
      {
        using CellShapeTag = vtkm::CellShapeTagQuad;
        f(CellShapeTag{}, std::forward<Args>(args)...);
      };
      break;
      default:
      {
        //EXPECT_TRUE(false); // 不希望运行到这里
      }
      break;
    }
  }
  template<typename PointVecType>
  VTKM_EXEC FloatDefault CellArea(vtkm::CellShapeTagTriangle&& shape,
                                  const PointVecType& points) const
  {
    using ShapeTag = vtkm::CellShapeTagTriangle;
    using QBase = Quadrature<ShapeTag, libMesh::Order::FIRST>;
    using FEBase = Lagrange<ShapeTag, libMesh::Order::FIRST>;
    FloatDefault area = 0;
    const auto& weights = QBase::weights();
    const auto& qpoints = QBase::points();
    for (int qp = 0; qp < QBase::NUM_POINTS; qp++)
    {
      auto local = qpoints[qp];
      auto jacobian = FEBase::jacobian(points, local);
      auto det = vtkm::MatrixDeterminant(jacobian);
      auto JxW = det * weights[qp];
      area += JxW;
    }
    return area;

    //auto directed_area = vtkm::Cross(points[0] - points[1], points[2] - points[1]);
    //return vtkm::Magnitude(directed_area) / 2.;
  }

  template<typename PointVecType>
  VTKM_EXEC FloatDefault CellArea(vtkm::CellShapeTagQuad&& tag, const PointVecType& pts) const
  {
    using ShapeTag = vtkm::CellShapeTagQuad;
    using QBase = Quadrature<ShapeTag, libMesh::Order::THIRD>;
    using FEBase = Lagrange<ShapeTag, libMesh::Order::FIRST>;
    FloatDefault area = 0;
    const auto& weights = QBase::weights();
    const auto& qpoints = QBase::points();
    for (int qp = 0; qp < QBase::NUM_POINTS; qp++)
    {
      auto local = qpoints[qp];
      auto jacobian = FEBase::jacobian(pts, local);
      auto det = vtkm::MatrixDeterminant(jacobian);
      auto JxW = det * weights[qp];
      area += JxW;
    }
    return area;

    //auto area = vtkm::Magnitude(vtkm::Cross(points[0] - points[1], points[2] - points[1]));
    //area += vtkm::Magnitude(vtkm::Cross(points[2] - points[3], points[0] - points[3]));
    //return area / 2.;
  }

  template<typename PointVecType, typename CellShapeTag>
  VTKM_EXEC void operator()(CellShapeTag& shape_generic,
                            const PointVecType& pts,
                            FloatDefault& area) const
  {
    auto cell_op = [&](auto&& tag)
    {
      using Tag = decltype(tag);
      area = CellArea(Tag{}, pts);
    };

    vtkm::CellShapeTagGeneric shape{ shape_generic.Id };
    OnCellOperator(shape, cell_op);
  }
};


void CalculateArea::Init() {}

void CalculateArea::solve()
{
  auto& mesh = static_cast<poisson::Mesh&>(get_mesh());
  mesh.prepare_for_use();
  auto data_set = mesh.GetDataSet();

  std::cout << "单元数：" << data_set.GetCellSet().GetNumberOfCells() << std::endl;
  std::cout << "节点数：" << data_set.GetCoordinateSystem().GetNumberOfValues() << std::endl;

  vtkm::cont::Timer timer;
  timer.Start();


  ArrayHandle<FloatDefault> area;
  vtkm::cont::Invoker{}(CalAreaWorklet{},
                        data_set.GetCellSet(),
                        data_set.GetCoordinateSystem().GetDataAsMultiplexer(),
                        area);

  auto toal_area = vtkm::cont::Algorithm::Reduce(area, 0.0);

  std::cout << "时间开销：" << timer.GetElapsedTime() << std::endl;
  std::cout << "total area：" << toal_area << std::endl;
}

